# TODO: Add comment
# 
# Author: roger
###############################################################################

library(quadprog)
library(tseries)

smi <- get.hist.quote("^ssmi",quote="AdjClose",compression="m")
smi <- smi[-dim(smi)[1],]

nesn <- get.hist.quote("nesn.vx",quote="AdjClose",compression="m")
nesn <- nesn[-dim(nesn)[1],]


x <- arima.sim(model=list(ar=c(0.5,-0.4,0.2),ma=c(0.9,0.1)),10000)
x <- as.numeric(x)
x <- as.numeric(diff(log(smi)))
x <- as.numeric(diff(log(nesn)))

x <- x[(length(x)-63):length(x)]
x <- x[(length(x)-(52*7)):length(x)]


arMax <- 12
n <- length(x)

usedSeq <- n - arMax

idxOrg <- 1:usedSeq
idxExpl <- t(t(matrix(rep((idxOrg),arMax),ncol=arMax)) + 1:arMax)
xr <- rev(x)
y <- xr[idxOrg]
X <- matrix(xr[idxExpl],ncol=arMax)

w <- rev(1:usedSeq)/sum(1:usedSeq)*usedSeq

Dmat <- t(X)%*%X
dvec <- t(X)%*%y

Dmat <- t(X*sqrt(w))%*%(X*sqrt(w))
dvec <- t(X)%*%(y*w)

Amat1 <- diag(arMax)
Amat2 <- diag(arMax)
diag(Amat2[2:arMax,1:(arMax-1)]) <- -1
Amat2 <- Amat2[,1:(arMax-1)]

bvec1 <- rep(0,arMax)
bvec2 <- rep(0,arMax-1)
		
Amat <- cbind(Amat1,Amat2)
bvec <- c(bvec1,bvec2)

#Amat <- cbind(Amat1)
#bvec <- c(bvec1)

#min(-d^T b + 1/2 b^T D b)
sol <- solve.QP(Dmat, dvec, Amat, bvec, meq=0) 

par(mfrow=c(2,1))
barplot(sol$solution)
barplot(sol$unconstrained.solution)

b <- sol$solution
buc <- sol$unconstrained.solution

fitted <- X%*%b
resid <- y-fitted
SS <- t(resid)%*%resid
sigma2 <- SS/(length(resid)-arMax)

fitteduc <- X%*%buc
residuc <- y-fitteduc
SSuc <- t(residuc)%*%residuc
sigma2uc <- SSuc/(length(residuc)-arMax)


fitteduc <- (X*sqrt(w))%*%buc
residuc <- y * sqrt(w) - fitteduc
SSuc <- t(residuc)%*%residuc
sigma2uc <- SSuc/(length(residuc)-arMax)


x0 <- y[1:arMax]

mu <- t(x0)%*%b
sigma <- sqrt(sigma2*(1+t(x0)%*%solve(Dmat)%*%x0))

muuc <- t(x0)%*%buc
sigmauc <- sqrt(sigma2uc*(1+t(x0)%*%solve(Dmat)%*%x0))

muuc + qt(c(0.025,0.5,0.975),length(residuc)-arMax)*sigmauc

print(mu)
print(muuc)
print(sigma)
print(sigmauc)
print(mu/sigma)
print(muuc/sigmauc)





